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We have studied the dynamics of a generahzed toric code based on qudits at finite temperature 
by finding the master equation coupling the code's degrees of freedom to a thermal bath. As a 
consequence, we find that for qutrits new types of anyons and thermal processes appear that are 
forbidden for qubits. These include creation, annihilation and diffusion throughout the system code. 
It is possible to solve the master equation in a short-time regime and find expressions for the decay 
rates as a function of the dimension d of the qudits. Although we provide an explicit proof that the 
system relaxes to the Gibbs state for arbitrary qudits, we also prove that above a certain crossover 
temperature, qutrits initial decay rate is smaller than the original case for qubits. Surprisingly this 
behavior only happens with qutrits and not with other qudits with d > 3. 
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I. INTRODUCTION 

It is known that the fragility of quantum states in the presence of interaction with an environment represents 
the main challenge for the large scale implementation of quantum information devices in quantum computation and 
communication. Quantum error correction is the theoretical method that was devised to protect a quantum memory or 
communication channel from external noise [1-8]. In these quantum error correction schemes, to improve the stability 
of quantum information processing, the logical qubits should be implemented in many-particle systems, typically TV 
physical spins per logical qubit. This is the quantum version of the classical method based on encoding information 
by repetition or redundancy of logical bits in terms of physical bits [9, 10]. The logical qubits should be stable objects 
with efficient methods of state preparation, measurements and application of gates. By efficiency we mean certain 
scaling behavior, e.g. the lifetime of a logical qubit should grow with N. 

In order to implement fault-tolerant methods for quantum information processing, we need to find a physical system 
with good enough properties to accomplish this protection from noisy environment and decoherence. One promising 
candidate are topological orders in strongly correlated systems. Here, the ground state is a degenerate manifold of 
states whose degeneracy depends on the topological properties of a certain lattice of qubits embedded into a surface 
with non-trivial topology [11]. Many-body interacting terms in a Hamiltonian are responsible for the existence of this 
topological degeneracy. The logical qubits are stored in global properties of the system represented by non-trivial 
homological cycles of the surface. In this topological codes, the property of locality in error detection and correction is 
of great importance both theoretically and for practical implementations [12-14]. It is also possible to generalize this 
topological codes for units of quantum information based on multilevel systems known as qudits, i.e. d-level systems, 
[15-18] and study its local stability [19]. An alternative scheme to manipulate topological quantum information is 
based not in the ground-state properties of the system but in its excitations [11]. These are non-abelian anyons 
that can implement universal gates for quantum information [20]. Yet, being within the framework of topological 
codes based on ground state properties, it is possible to formulate new surface codes known as topological color 
codes (TCCs) [21] such that they have enhanced quantum computational capabilities while preserving its nice locality 
properties [21-23]. TCCs in two-dimensional surfaces allows for the implementation of quantum gates in the whole 
Clifford group. This makes possible: quantum teleportation, distillation of entanglement and dense coding in a fully 
topological scenario. Moreover, with TCCs in 3D spatial manifolds is possible to implement the quantum gate tt/S 
thereby allowing for universal quantum computation [24, 25]. Very nice applications of topological surface codes have 
been found in other areas [26, 27]. 

Acting externally on topological codes, in order to cure the system from external noise and decoherence, produces 
benefits from the locality properties of these codes. Namely, a very important figure of merit is the error threshold 
of the topological code, i.e. the critical value of the external noise below which it is possible to perform quantum 
operations with arbitrary accuracy and time. For toric codes with qubits, the error threshold is very good, about 11% 
[12]. This value is obtained by mapping the process of error correction to a classical Ising model on a 2D lattice with 
random bonds. Interestingly enough, this type of mapping can be made more general and applied to TCCs yielding 
the same error threshold [28] while maintaining its enhanced quantum capabilities [29, 30]. These results have been 
confirmed using different types of computation methods [31-35]. It is also possible to carry out certain computations 
by changing the code geometry over time, something called 'code deformation' [12, 36, 37] that allows us to perform 
quantum computation in a different way. A more general type of codes can be constructed with quantum lattice 
gauge theories based on quantum link models [38]. 
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In this paper we adopt a different approach than external protection of topological codes. Hence, instead of 
performing active error correction, we just rely on the robustness of a Hamiltonian that has a gap above the ground 
state manifold where the quantum information is stored. Thus, we leave the system to interact with the surrounding 
environment and study the fate of the topological order under these circumstances. This source of noise is inescapable: 
the microscopic interactions of the physical spins with thermal particles or excitations of the local environment. The 
analogous situation for classical information processing is well-understood, but the existence of a similar mechanism 
for quantum information is still an open problem. The quantum theory of open systems provides a natural framework 
for studying stability in the presence of thermal noise. The particularly simple properties of Kitaev's model allow 
us to apply the Davies' theory, namely the dynamics of a quantum system weakly interacting with a heat bath in 
the Born-Markov approximation [39-46]. There have also been some related studies regarding thermal effects on 
Adiabatic Quantum Computation [47, 48]. 

The first indication that the toric code for qubits in 2D spatial dimensions is unstable against thermal noise 
was shown in [12]. Further analytical and quantitative arguments of thermal instability were given by [49]. Later, a 
rigorous proof of this fact has been estabhshed using the theory of quantum open systems [50, 51]. Subsequently, other 
investigations have been conducted for abelian models, non-abelian models, TCCs [52-56] etc. Remarkably enough, 
while with qubits in 2D lattice models the topological protection is lost under the action of thermal fluctuations [57] , 
however it is possible to set up a full-fledged topological quantum computation using certain types of topological color 
codes in higher dimensional lattices [58]. Under these conditions, it is possible to prove that self-correcting quantum 
computation including state preparation, quantum gates and measurement can be carried out in the presence of the 
disturbing thermal noise. Additionally note that thermal noise does not always turn out detrimental in quantum 
information, even for systems without topological order [59, 60]. 

In this work we extend those results regarding the thermal effects on generalized toric codes constructed out of 
qudits. We hereby summarize briefly some of our main results: 

i/ We formulate the dynamics of a generalized toric code based on qudits at finite temperature. To this end, we find 
the master equation coupling the qudits of the system code to a thermal bath. 

ii/ We study and classify the different types of thermal processes that may occur when the anyonic excitations are 
created, annihilated or diffused throughout the system. In particular, we find that for qutrits new types of anyons 
and thermal processes appear that are forbidden for qubits. 

iii/ The master equation is too involved so as to yield an explicit expression for the decay rate of the topological 
order initially present in the code. However, in a short-time regime it is possible to solve it and find expressions for 
the decay rates as a function of the dimension d of the qudits. Interestingly enough, we find that the decay rate for 
qutrits presents a crossover temperature Tc that is absent for any other qudits. 

iv/ We can give an explicit proof that for times long enough, the non-local order parameter representing the topological 
order in the system decays to zero. 

This paper is organized as follows: In Sect. II we review the formulation of the master equation of the 2D Kitaev 
code for qubits in order to establish the notation and the necessary tools to study thermal effects in more general 
toric codes. We also introduce a non-local order parameter and study the fate of topological orders for two different 
regimes: short time-regime and long-time regime. In Sect. Ill we find the master equation describing topological 
qutrits coupled to a thermal bath. This allows us to see new energy processes for the anyonic excitations that are not 
present when the toric code is made up of qubits. Likewise, the short-time regime has a different behavior that can be 
seen in the initial decay rate of the topological order. In particular, we can define a crossover temperature for qutrits 
where the decay rate is better than with other qudits. Sect. IV is devoted to conclusions. We refer to Appendix A 
for evolution of the order parameter for qutrits and Appendix B for a proof of the irreducibility of the computational 
representation of the (i-Pauli group needed to study the master equation in the long-time regime. 

II. THERMAL STABILITY OF THE KITAEV 2-D MODEL 

We shall not dwell upon all details of Kitaev's toric code [4], however we will give the basic ideas to understand 
how to apply a thermal stability analysis to it, as well as to establish the notation and methods. We will consider a 
k X k square lattice embedded in a 2— torus. Let us attach a qubit, like a spin 1/2, to each edge of the lattice. So we 
have N = 2k'^ qubits. For each vertex s and each face we denote the stabilizer operators of the following form: 

A, := n ■= n (1) 

jGstar(s) jGboundary(p) 

where Xj and Zj are the Pauli Matrices applied to the qubit on site j. As and Bp commute among each other for they 
have either or 2 common edges. They are also Hermitian and have eigenvalues 1 and —1 (see figure 1). Therefore, 
they constitute an abelian subgroup of the Pauli group of n qubits that is a stabilizer group. 
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FIG. 1: Square lattice on the torus. The yellow points represent qubits. 



Let 1-L be the Hilbert space of all n = qubits and define the topological quantum code or protected subspace 
C C H as follows: 

C = [\^)en: ^|^) = |^), 5p|^) = |^)forall5,p}. (2) 



This construction defines a quantum code called the toric code. The operators A^, Bp are the stabilizer operators 
of this code, i.e. operators that leave trivially invariant the code space. As we want to analyze the physical properties 
of this code, in particular the thermal properties of the topological order, it is convenient to define its associated 
Hamiltonian in the form: 

S p 

Complete diagonalization of this Hamiltonian is possible since operators Ag, Bp commute. In particular, the ground 
state coincides with the protected subspace of the code C; it is 4-fold degenerate (see figure 2). All excited states are 
separated by an energy gap AE > 4. This is due to the fact that the difference between the eigenvalues of As {Bp) 
is equal to 2. Excitations come in pairs since they correspond to violations of the plaquette and/or vertex stabilizer 
operators and these must comply with the overall constraints Hs = ^ Hp = 1- Thus, excitations are 
represented as open strings in the direct or the dual lattice of the original square lattice. 

Energy 

yv 2'' - 4 

Excited states 




4-degerate ground state (CODE SPACE) 

FIG. 2: Schematical spectrum of the Toric Code Hamiltonian. The ground state is the code space C where we codify our 
information. 

An essential feature of this Hamiltonian is its locality in terms of four-body interactions, very useful for practical 
purposes. Another key property is that this Hamiltonian model is gapped, which led to the initial expectation that all 
type of "errors" , i.e. noise-induced excitations will be removed automatically by some relaxation processes. Of course, 
this requires cooling, i.e. some coupling to a thermal bath with low temperature (in addition to the Hamiltonian (3)) 
as we shall describe later on. It can be shown that this Hamiltonian is robust under local quantum perturbations at 
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zero temperature [57]: there would be a level splitting which will vanish as exp(— a/c), where k is the length of the 
lattice [4]. 

Due to this unavoidable coupling to a thermal bath, our system is subjected to thermal errors as well. These can 
be seen as violations on the plaquette and vertex conditions: = |^), = |^)- Moreover, As and Bp are 

unitary, and also Hermitian in the case of qubits. Therefore, violations on the plaquette and/or vertex condition mean 
are given by 

A,\^) = -I*), Bp\^) = -I*), (4) 

for a certain number of sites s and/or plaquettes p. 

These violations cost energy to our system, thereby becoming excitations. And as long as they always come in pairs 
(to satisfy the condition Yl^ Ag = I and Yl^ Bp = 1), they can be seen (pictorially) as string operators with plaquette 
or vertex violations at the ends. 

Errors on the system can be expressed in terms of operators a^, a^, or products among themselves. These operators 
act on each edge j where the physical qubits are placed. We use the notation for a Pauli operator of type X when 
it is referred to an error, i.e., a bump operator acting due to the coupling to the thermal bath. Similarly with . 
Namely, it is just a matter of notation to distinguish when we have an operator that defines our stabilizer operators 
in As, Bp and when we have an error acting on the system. To see what is the effect that they produce, we will see 
how the ground state changes by applying these cr^'^. We will see that this corresponds to creation, annihilation and 
movement of a pair of excitations, that from now on we shall refer to them as anyons. They are called anyons since 
their wave function picks up a different phase than fermions or bosons when we exchange the end-particles of string 
operators of x-type with z-type. According to this notation, when we apply a bump operator from the thermal bath, 
it will act on the ground state of the system as follows: 

(5) 

where |^) is the ground state of the system where our information is encoded. This means that the physical qubit 
at the edge j has been bumped. The energy cost will be AE = 4 in energy units of the system corresponding to the 
definition of H^^^ . 

As a first step, one is interested in designing a stable quantum memory, i.e. a A^— particle system which can support 
at least a single encoded logical qubit for a long time, preferably with this time growing exponentially with N. This 
is the notion of stability that we shall refer to from now on. In the paper [50] by Alicki et al. they provide a rigorous 
method to prove thermal instability of the 2D Kitaev model and obtain a master equation that describes the dynamics 
of the system weakly coupled to a thermal environment. We will study the problem of thermal instability within the 
framework of topological orders obtaining complementary and interesting results. 

A. Davies' Formalism 

Let us consider a small and finite system, that is coupled to one or more heat baths at the same inverse temperature 
P = (/cbT)"-^ leading to the total Hamiltonian 

H = H'^' + H^^'^ + V with y = ^ 5'a /a. (6) 

Here H^^^ represents the Hamiltonian of the system, where the quantum information is encoded and which we want 
to protect from the external thermal noise, i^^ath ^j^^ bath Hamiltonian, i.e., it describes the internal dynamics of 
the bath which is out of our control. Finally V represents the coupling between the system and the thermal bath. 
Sa and fa are operators which act on the system and bath respectively. Both the coupling operators So, and are 
assumed to be Hermitian (without loss of generality [41]). 

In the weak coupling regime that we shall assume throughout this work, the Fourier transform of the auto- 
correlation function of fa plays an important role as it describes the rate at which the coupling is able to transfer 
energy between the bath and the system [39-42]. Often a minimal coupling to the bath is chosen, minimal in the 
sense that the interaction part of the Hamiltonian is as simple as possible but still addresses all energy levels of the 
system Hamiltonian in order to have an ergodic reduced dynamics. This last condition is ensured if [41, 43-46] 

{S„, H^y^y = C1, (7) 
i.e. no system operator apart from those proportional to the identity commutes with all the and H^^^. 
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The weak coupling limit results [39-42] in a Markovian evolution for the system given in Heisenberg picture by the 
master equation 

dX 

— =g{X):=iS{X) + C{X). (8) 

The generator of the evolution g{X) is a sum of two terms, the first is a usual Liouville-von Neumann term as in the 
quantum mechanics of closed systems, while the second is a particular type of Kossakowski-Lindblad generator: 

S{X) = [H''y\X] (9) 

a uj>0 
a uj>0 

+ e-^^Sa{io) [X, (5„H)^] +e-^^ [S^{io) , X] (5„H)^}. (11) 
Here the Sa{ui) are the Fourier components of Sa as it evolves under the system Hamiltonian 

^^'"'-S^e--''"'''=Y^S^iu:)e---\ (12) 



e 



where the cj's are Bohr frequencies of the system Hamiltonian {huj = Ei — £^2, for two energy levels Ei and £^2)- 

In addition, the temperature of the environment appears in (11) through /3, and this generator is the so-called 
Davies generator [39], or Born-Markov generator in the quantum optics literature. 



B. Master Equation for 2-D Kitaev Model with Qubits 



Given the simplicity of the Kitaev's model, we can apply the Davies' theory for studying its stability in the presence 
of thermal noise. This is pictorially represented in figure 3. 



HEAT BATH 




FIG. 3: Toric Code coupled to a heat bath. Outgoing arrows in the upper part of the figure mean information fiowing from 
system-to-bath, and ingoing arrows in the lower part mean information fiowing from bath-to-system. 

The interaction Hamiltonian is assumed to be local and associated with and errors: 

y = Y.^i^f!+^j®f!^ (13) 

j 

where fj and /J are associated with two different baths. Thus, first of all, we need to compute the Fourier transform 
of the system operators e^^^'^^Vje"^^^'^^'' and e^^^"^^^ (j| e-^tH^^'' order to define the dynamical operators of the 
system. Here H^^^ := Hq = - Yls - Ep Bp, with [As, Bp] = 0, [A^, cr|] = and [Bp, aj] = 0. Thus, stabilizers As 
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only play a role in the Fourier transform of <jj and Bp only in <j|. By computing this Fourier transform, we obtain 
the dynamical operators of the system due to the coupling to the thermal bath. With A = 4 denoting the gap of the 
Toric Code Hamiltonian, then the expression of these operators Sa{oj) that appear in equation (11) are [50] : 



1. Operators associated with errors: 



SJ{0) 


■= 


'£ pO 


5J(A) 


■■= h 




SJi-A) 







with i?^ := ^(1 — BpBp') and := ^{1 ^ Bp){l ^ Bp') being orthogonal projectors. 
2. Operators associated with cr| errors: 



S]{0) 








:= aj 




S^i-A) 


:=a] 





and the projectors: := ^(1 - A^As^) and P^ := i(l ^ As){l T As'). 



These operators have a nice interpretation in terms of anyonic properties of the system: 

1. ctj{b^j) creates a pair of anyons of z-type(x-type) on the lattice at position j. See figure 4. 

2. ajihj) annihilates a pair of anyons of z-type(x-type) on the lattice at position j. See figure 5. 

3. oPjip^) moves a pair of anyons of type (x- type) on the lattice. See figure 6. 




FIG. 4: Dynamics induced by the heat bath on the Toric Code. Creation of a new pair of anyons. Energy increases — 4. 




FIG. 5: Dynamics induced by the heat bath on the Toric Code. Annihilation of a pair of anyons. Energy goes down by /\E — 4. 
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FIG. 6: Dynamics induced by the heat bath on the Toric Code. Pure decoherence by moving an anyon with no energy change. 



Thus, the dissipator of the master equation C{X) for the system is: 



2 



+ h,Xh]}-\Rmhl[h%X]i 
C\X) = -R{4){{-a]ajX - Xa]aj + a]Xaj) + e-^^{-aja]X - Xaja] 



2 

+ ajXa]}-^R{0)[al[alX]]. 

Where i?(4) and R{0) are the exchange rate between the system and the bath associated to each Bohr frequency, 
namely uj = 0, 4, assuming units of J = 1. 



C. Topological Order 

We shall study the evolution of the expectation value (GS| Xc |GS) as a simple order parameter, where Xc is the 
tensor product of Pauli operators along one non-contractible loop on the surface of the torus and |GS) denotes a 
generic ground state of the system Hamiltonian. This ground state is a superposition of the degenerate states in the 
ground state manifold of i^^^®, namely C. This gives us a sufficient measure of the topological order of the system [49]. 
If this quantity falls to zero during the time evolution for every element of C, there is not a global and self-protected 
way to encode quantum information. The evolution of the operator Xc is given by equation (47), 

In order to simplify the computation, we remove the free evolution by performing the transformation 

X^{t)=e--^"''''XS)e^"'''\ (17) 
Since the dissipator is invariant under this transformation, we obtain 

dXe{t) 



dt 



C[X,{t)]. (18) 



Interestingly, for the expectation value we obtain (GS| Xc{t) |GS) = (GS| Xc{t) |GS), as |GS) is an eigenstate of H^^^ . 
Taking into account expressions (14) and (15) , the action of the dissipators on Xc can be simplified to 

C-.{Xc) = -^E^(^) {[RidRi,Xe]]+e-^^[B:^_,[iy_,X,]])+R{0)[Ri,[Ri,X,]], (19) 
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and 

CMc) = ^i?(A)[Pi<7|Xe<7|Pi - PiXc + e-^^(Pi<7|Xe(7|Pi - PLX,)] 

3 

+ P{0)[i^a|Xea|i^-P^XJ, 

where we have used the fact that [P^ -^c] = for every j, as these projectors are only functions of vertex operators. 

However, the same assertion is not true for i?^ q in general. If j c, i.e. j does not belong to the path where Xc is 
acting on, every element commutes with each other and their contribution is zero. On the other hand, if j G c, as 
a^a^aj = — crj, the string operator yields a^XcCTj = —Xc- Therefore, simplifying we obtain 

= -^\c\X,{R{A)[Pi+e-^''Pi]+R{0)Pi}, (20) 
where |c| is the number of points in the path c. 



D. Short-Time Regime 



The solution to the master equation (18) is formally written as Xc{t) = e^^Xc. However, this expression is too 
involved to be computed analytically except for short and long times to be specified hereby. In the first case, at lowest 
order we have 



Xe(t)^(l + t£)Xe. 

The evolution of (GS| Xc(t) |GS) is given by 

(XeW) ^ [1 - 2t\c\R{A)e-^^]{Xcm- 
To arrive at this equation, we have used the fact that for all j: 

n,o |GS) = 0, 
Pi\GS) = |GS), 

P^^oIGS) = 0, 
Pi|GS) = ICS). 

Thus, the contribution of is zero: 

(GS|£,(X,)|GS) = -lj2R{A)e-'^^ {GS\[RL,[Ri,X,]]\GS) 



1 ^P(A)e-^'' ((GS| |GS) - (GS| [Pi,Xe] |GS)) = 0. 



(21) 



(22) 



(23) 



Whereas for Cz, we have 



£,(Xe) = --|c|P(A)e-^^ (GS| X, |GS) . 



Finally, as (GS| Xc{t) |GS) = (GS| Xc{t) |GS), the desired equation valid at short times is 

A 



1 - =i|c|P(A)e-^'^ 



(24) 



(25) 



with A = 4. 

It is important to remark that P(0) does not appear in the initial decay rate, as long as short times is concerned. 
The diffusion of anyons is a second-order process in time as it requires first the creation of a pair of anyons with P( A) , 
and later the free diffusion with R{0). 
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E. Long-Time Regime 

On the other hand, in order to analyse the thermal properties for long times, we write the Davies generator in the 
Schrodinger picture through the relation Tt[C^ {p)X] = Tr[p>C(X)] for any X and p. It is a well-known result [39-42] 
that the Gibbs state is a stationary state for , 

C\P0) = O, (26) 

where = e~^^^^^/Z, P is the same to the inverse temperature as the surrounding bath, and Z := Tr(e~^^^^^) 
is the system partition function . To guarantee that any initial state of the system relaxes to we can resort to 
condition (7). In our case this follows from the Schur's lemma as S'q, = crj, cr| and {1, cr^, cr^, a^a^} form an irreducible 
representation of the Pauli group. 

Thus (GS| Xc{t) I OS) ^ Ty[XcP/3] for large t, and we have Tr[XcP/3] = 0. This is simply due to the fact that p^ is 
diagonal in any of the possible eigenbasis of H^^^, and it is not difficult to choose one such that Xc vanish on diagonal 
elements, 

i 

for some eigenbasis oi Hq, the Kitaev's Hamiltonian. 

In conclusion whatever the initial value of the order parameter (Xc(0)) is, it decays to zero during the time evolution 
of the system provided that the temperature is finite. The decay rate at short times is equal to ^\c\R{A)e~^^ . Note 
the detrimental effect of the factor |c|: the larger size of the system, the larger the decay rate. In order to keep the 
order parameter above certain finite value such that {Xc{0)) ^ 0, this decay rate must decrease, which is not the case 
when increasing the system size. 



III. KITAEV 2-D MODEL FOR QUDITS 

In this section we consider again a 2-D toric code, but instead of assuming that we have a two-level system on 
each site, we will consider that particles arranged on the torus, have d accessible levels. We will first derive a general 
theory for qudits and then consider the case d = 3 (qutrits). A qutrit can be represented for instance, as a particle of 
spin 1 or a 3-level system in an atom, etc. 

This problem is very interesting since qutrits have certain advantages with respect to qubits, namely: 

1. Larger capacity of information storage. 

2. Quantum channels are more robust for qutrits. For example: Bell inequalities are proved with more accurate 
bounds. This is relevant for Quantum Key distribution. 

3. Entanglement quantum destination is more efficient with qutrits than with qubits [61]. 

4. Qutrit logic gates [62] are also capable of providing universal quantum computation, i.e., the necessary compu- 
tational power to construct all possible logic gates [8]. 

To build a system like that, we will try to choose the Hamiltonian and the operators acting on the system in the 
same way as before. Previously, for two-level systems, we have considered the Pauli matrix algebra to be the basis of 
operators in our system. Now, we have to use a proper generalization for dimension d. As iXZ = Y gives the second 
Pauli matrix, it is enough to consider X and Z in this generalization to quantum states with d multilevels. However, 
the generalization of Pauli matrices to dimension d is not unique [63]. Thus, we shall select the most important 
properties of Pauli matrices of dimension 2 for our purpose of quantum error correction. 

In d = 2 we defined a basis: |0) , |1) in the Fock space of each particle. They are defined as the eingenstates of the 
Z Pauli matrix. And the X Pauli matrix takes |0) to |1) and viceversa. 



Z\0) = 
X\0} = 



n 




1 y ' 




+ |0), 


Z\l) 


+ |i>, 


X\l} 



= -ll) 

= |0), 



The key important properties of these matrices for doing error correction are: 
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• They satisfy a cyclic condition (i.e. applying twice Z or X Pauli matrices is the identity) i.e. they are unitary. 

• They anticommute, that means XZ = —ZX. 

Those are the properties that are generalized to the d— dimensional case. Hermiticity is not taken into account as 
a basic ingredient, as we can always add the Hermitian conjugate obtaining an Hermitian operator, e.g. Z = Z ^ Z\ 
then Z is Hermitian. Now we consider a basis for the particle Fock space: |0) , |1) \d— 1) which will be the 
eingenvectors of the generalized Z matrix with a certain eingenvalue. We define X as the operator which takes the 
state |0) to |1) then |1) to |2) and so on. We will also ask for a cyclic condition as in the previous case: 

= 1, = 1. (28) 

All these requirements can be cast on to the following defining relations: 

Z |0) = + |0) , Z |1) = |1) Z |2) = |2) , ... , Z |d - 1) = cj^-^ \d-\)\ 
X|0) = X|l) = |2), ... ,X|d-l) = |0). 

(29) 

Looking at equation (29) we can deduce the meaning of operators X and Z. X is the displacement operator in the 
computational basis (i.e. in the Fock space basis of the physical qudits). Z is the dual operator of X under a discrete 
Fourier transform. In other words, Z is diagonal in the computational basis and its eingenvalues are the weights of 
the Fourier transform. Thus, X plays the role of the displacement operator and Z is the dual operator on a system 
with discrete states of qudits [8]. 

Due to the cyclic condition (28) of Z (Z^ = 1) we have the relation cj^ = 1 where in general is a complex number. 
This implies that cj is a primitive (i— root of the unity, 

cj = e^^. (30) 

Additionally, we can easily verify that ZX — ujXZ^ as follows from equation (29). 

We have already the algebra of operators that we are going to use in order to built the stabilizer operators on this 
qudit toric code. The problem is that if we construct the vertex and plaquette operators as before, namely, 

A, = n = n (31) 

jGstar(s) jG boundary (p) 

then [As, Bp] ^ for all s and p. They commute with each other provided that they do not share any common edge, 
but that is not the case if they share two. This happens because in this case the operators X and Z are no longer 
Hermitian. 
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FIG. 7: New lattice for qudits showing vertex As and plaquette Bp operators: orientation of the lattice is necessary. 
As shown in figure 7, we have 

[As, Bp] = [X1X2X3X4, Z3Z4Z5Z6] = (1 - cj^)AsBp, (32) 

which does no vanish for for dimension d > 2, (1 — cj^) 7^ 0. The case of d = 2 is a very special case with uj = —1 
and therefore (1 — cj^) = 0. This happens because for d = 2, X and Z are Hermitian operators. We need to think 
of another way to define our operators to have the same commutation rules as before, and this leads to define an 
orientation on the lattice. This is shown in figure 7. Defining an orientation on the lattice is a direct consequence of 
the non-Hermiticity of operators X and Z. 
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Using the orientation of the lattice, we define the stabilizer operators in the following way. To build the vertex 
operators ^4^ we assign an operator X or depending on the arrows of the edges of the lattice. If an arrow is 
pointing towards the vertex j, we will use XJ^ to build A^, and if the arrow is pointing out another vertex we use 

Xj. For plaquette operators Bp^ is taken if the arrow is pointing clockwise and Z^^ for anti-clockwise as shown 
in figure 7. To see now that we obtain the correct commutation rule, we look again at figure 7 and check, 

[As, Bp] = [Xf ^X2X3-^X4, Zg-^ZgZsZg-^] = (1 - ujuj-^)AsBp = (33) 

Then, the Hamiltonian could be written as follows: 

S p 

Although, according to the definition of As and Bp, this operator is unitary, it is important to note that the operators 
As and Bp are not Hermitian any more, so i^^aux is no longer Hermitian. However we may redefine the Hamiltonian 
in the following way: 

H^y^:=^{H,^^ + HlJ, (35) 

where H^^^ is Hermitian now. The effect that i?j[ux ^^le system is a redefinition of the orientation on the lattice. 

So we have a superposition of a lattice orientated in the way of figure 7 (arrows up and right) and another with 
arrows down and left. Nevertheless, one can always think in terms of i^aux for the pictorial image and then use H^^^ 
to compute energies and derive equations. 



A. Anyon Model 

The theory developed above was done for the general case of qudits. From now on and to be concrete concerning 
thermal effects, we will focus on the case where d = 3 (qutrits). Later on we will be able to extract conclusions for 
qudits as well. There are still many important aspects to be studied about this model and its coupling to a thermal 
bath. We need to compute the energy gap of the Hamiltonian, i.e. the energy difference between the ground state 
where the code lies and the excited states which represent the errors. It is also important to calculate the anyon 
statistics, as long as they are associated with the excitations of a topological system with qutrits. 

In d = 3 the phase factors are uj = e^^ , cj^ = e^^ , cj^ = 1. We will see for this particular case, how excitations can 
be created, moved and annihilated. This will give us the properties of the anyon model which is going to be associated 
with the group Z3. 

As we did before, we use a notation in which <j| = Zj and crj = Xj, except that we use the symbol a to denote 
errors acting on the system, i.e., bump operators acting because of the coupling to the thermal bath, whereas we shall 
use X, Z for the Hamiltonian interactions defined by the vertex and plaquette operators of H^^^. 

Errors on the system can be expressed in terms of operators a^, or products containing them, and acting on 
each edge j where the qutrits are placed. And the same goes for . To see what is the effect of these errors on the 
system, we will see how the ground state changes by applying cr^'^. We will see that this corresponds to processes in 
which anyons are created, annihilated or moved throughout the torus. 

Let us see what happens when we bump a qutrit in a position j from the outside and then act with the Hamiltonian 

-f^aux5 

Note that every operator of the Hamiltonian commutes with this cr| except two As operators which share a leg with 
this qubit j. But, contrary to the case of d = 2 there is an orientation defined on the lattice. So, for instance, if an 
error (cr|) occurs in a certain vertical edge, one of these As (the one below) is defined with an Xj, thus: 

A,a] \^P) = uj-'a'jA, |V) = wV| , (36) 

but the As' above the edge is defined with X~^, then: 

A,,a] IV) = oja^As' \^} = toa] \^} (37) 

Hence, we have two violations of the vertex condition, one with charge uj and the other with uj'^. This is one of the 
two types of anyons that we will have in this system, and we shall denote it as an u;^ — uo anyon. It is important to 
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point out that these are only labels to classify the excitations based on the violations of the operator As (and Bp). 
In principle we could classify anyons based on the violation of stabilizers Aj^ (and B~^) that appears in H^^^. It is 
just a matter of labeling, the physics is the same. 

Now we can act with (j| again and obtain the other anyon type called uj — cj^. Actually they could be considered 
as the same anyon type as before but with opposite orientation. However, it is convenient to define them as two types 
of anyons as they will have different braiding properties. Moving anyons of the same type around each other will be 
different from the case of having anyons of different type. Likewise, it will be necessary to have anyons of different 
types in order to have fusion of anyons without annihilation. We shall explain this in the next subsection in more 
detail. 

Note that by acting twice with cr| is equivalent to act with (cr|)~^. Thus, although every error can be expressed in 
terms of X and Z operators, it will be useful to think sometimes as if we act either with X, Z or Z~^. All these 
arguments are exactly the same in the case of Bp operators and errors. Therefore, we have 4 types of anyons, 2 of 
plaquette type and 2 of vertex type. 
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FIG. 8: Anyons of type x (red) on the direct lattice. Anyons of type x (green) on the dual lattice. 

Let us study now the braiding of the anyons. We will consider two chains of different type: plaquette anyon and 
vertex anyon (as in figure 8). In this case we get something remarkably different from the d = 2 case. Now it is not 
the same to let one anyon still and move the other around it than do it the other way around. Thus, let us move 
particles around each other. For example, let us move an x-type particle around a z-type particle (see figure 9). 
Then, 

l^initial) = S'{t)\r{q)) , l^final) = S^c) {t) {q)) = | ^initial) , 

because S^{c) and S^{t) cross each other just on one qutrit satisfying the relation 

XZ = uj'^ZX 

and S'^(c)|'0^(g)) = |'0^(g)). We see that the global wave function, i.e. the state of the entire system, acquires the 
phase factor cj^. Nonetheless, if the operation is the opposite, that is, if we move a z-type particle around a x-type 
particle then: 

|*i„itial) = S%q)\r(t)) , |*fi„al> = S'{c)S%q)\r{t)) = C^|*initial>, 

since S^{q) and S^{c) cross each other just on one qutrit again satisfying the relation: 

ZX = uXZ (38) 

and S^{c)\ilj^{t)) = \ip^{t)). We see that the global wave function acquires now the phase factor uj. 

Therefore, we arrive at a very important novelty for qutrits that is different from when we dealt with qubits 
regarding two aspects: 

1. The phase that the anyon picks up is different from —1. 

2. The phase depends on the orientation in which the braiding close path is traversed. 
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FIG. 9: Anyons of type Z (red) on the direct lattice attached to a string t. Anyons of type x (green) on the dual lattice 
associated to a string q. The x-type particle moves around a z-type particle on a closed string c. 



B. New Anyon Energy Processes 

First of all, let us look at the gap of the Hamiltonian. We will reach our first excited state by applying a or 
operator to the ground state. Let us see which is the energy difference between the ground state and the first excited 
state. Remember that 211^^^ = i^aux + ^lux- We denote P and S the number of plaquette and vertex operators 
respectively, with P -\- S = N the number of qutrits in the lattice, and {l^l'} are the adjacent vertices of the site of a 
qutrit j: 

= i{- E ^ - E Bp + h.c} 1^) = -{P + S) |^> (39) 

S p 

Ha^lij) = ^-{-J2a, -J2Bp + h.c}a]\ij) = -{P + S-2)a]\ij)-^{Aia]\^p) + 

S p 

- A,a] l^P) + A\a] |V) - A\,a] m = -{P + S - 2)a| - wV| \^P) - Loa]4 ^ = 
= -{P + S-2 + L0 + LO^)a] \i^) = -{P + S-2 + 2 cos —)a] |V') = 
= -{P + S-3)a]\i>). 

Thus, the energy difference is 



AE = 3. 



The action of produces the same energy increment but we have to do the commutation with the operators Bp. 
This calculation can be easily extended to the case of qudits with arbitrary (i, obtaining the gap equation 

= Arf = 2 ^1 - cos . (40) 

Note that there is a reduction of the energy gap for d = 3 in comparison with the case of qubits, where it was 4. 
It is also important to point out that if we act again on the same bond of the lattice with (<j^)~^, there would be an 
energy reduction of the same amount of energy. Moreover, if at the endpoint of an anyon cu — cj^ we act with 
we obtain the same pair of anyons again, and same energy, but longer (see figure 10.2). In this process the energy is 
preserved AE = 0. This means that there is no energy exchange between the thermal bath and the system. We can 
understand the process as a diffusion of the anyon with no energy cost. In analogy to the case d = 2, this is what is 
called moving an anyon. It is also important to remark that still for qutrits, all process that involve moving a simple 
pair of anyons have no energy cost. 

Until here, there is a complete analogy with the case of d = 2. But we are going to see now a process that only 
occurs in d > 2. Imagine that there have been two excitations on the system, and two anyons of opposite orientation 
have been created. Moreover, they are separated by just one vertex operator. The situation is plotted in figure 10.1. 
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FIG. 10: 1) Fusion of anyons (ending tied, not annihilated). 2) Movement of an anyon. We plot just one dimension as long as 
the rest of the lattice is irrelevant, i.e, the process is the same everywhere. 

Imagine that we act now with a on the bond, which is error free, that links the anyons uJ^ — uj and uj — uP' 
(opposite orientation). Let us analyze the energy process. 

iiW) = -^(-5^^. - ^i?p)|^')+/i-c. = -(P + 5-6)|^), (41) 

S p 

Ha^jW) = -\{-Y.As-Y,Bp)a'^W)+h.c=-\{P + S-uj-uj-LO^LO^)a^W)- 

S p 

- ^{P^S-uo^-uj^- ujuj)ct' \^') = -(P + 5 - 6 + ^)ct' \^') , 
so, the energy difference is 

AE = -3/2. 

What has occurred is that two anyons have been tied together, but not annihilated. This process lowers the energy 
of the system in a smaller amount than the process of annihilation. If in this situation we would act with a on 
the point where the two pair of anyons are tied together, the two anyons would split apart, and this process would 
cost energy AE = 3/2. This could be analyzed exactly the same way with errors and Bp operators. 

It is remarkable that this phenomenon cannot happen in (i = 2, as in (i = 2 the product cucu = ( — 1)(— 1) = 1. 
Therefore, d = 3 is the first non trivial case to have processes like these in a toric code with qudits. 



C. Master Equation for Topological Qutrits 

As we have seen, all these processes are generated by the action of operators a^, (cr^)^ and a^, (cr^)^; as in this 
case, the square of the Pauli operators are their Hermit ian conjugate. Nevertheless, the energy exchange depends on 
the situation of the system when we bump it with the thermal bath from outside. Before writing the master equation 
that describes the dynamics of the system, it will be useful to distinguish between these situations by local projectors. 
The answer to the question whether this is possible or not in this case is not trivial. However, we show that it is 
possible to classify into groups of processes that have the same energy gain from the bath. Furthermore, they could 
be distinguished by certain projection operators that only involve two adjacent vertex or plaquette operators. 
We arrive at the following classification: 
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(42) 

In this table we have represented ah combinations of two adjacent topological charges. In the first column: we 
depict a representation of the different types of anyons, with two topological charges attached at their ends and linked 
by a dash. Correspondingly, all these anyons have an intrinsic orientation. At the left side of the dash there is the 
eigenvalue of the operator As and at the right side, the eigenvalue of the adjacent operator A^^. A physical qutrit j 
would be in the middle of the dash (see an example at figure 11). In the second column: we write the projector that 
gives 1 for that situation and for the others. 




FIG. 11: 1) Initial state 1 — 1 =^ Final state cj^ — u 2) Initial state cj^ — u =^ Final state cj — uo^. This is an example 
of what happens to the topological charges when there is a bump from the thermal bath outside. The first one the energy gain 
is = 3. The second one A^; = 0. 

Here we have defined the following operators in order to simplify the notation: 

4'2o,+i,-i(s,s') := 

AA{s,s') := A^^Us')-A^^Us), 

where 5 and 5' are the two vertex surrounding the qutrit j. The index a takes values on the exponent of the phases uj 
that appear from the braiding processes. These projectors tell us which are the charges of the system that surround a 
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certain qutrit. That is why they are local projectors. Moreover it is easy to verify that they form a set of orthogonal 
projectors: 

a. 

pj _ pit 

{pif = Pi- 

As we have already explained, we classify the situation of the system in terms of the charges according to the 
eingenvalues of the operators Ag associated with the part of the Hamiltonian i^aux- One could do the same thing 
for ^7^, but the situation of the system will be the same independently of the label we assign to them. So these 
projectors can discriminate perfectly between eigenstates of the Hamiltonian H^^^ . 

Now, given a certain state of the system l^/;'), by applying these projectors we can figure out which situation we 
have. This means that if an operator or (or their Hermitian conjugate) is going to act on our system, we will 
know which energy process is bound to happen. Based on this, and studying the different situations that we can 
encounter, one can define a set of operators that tells us whether an anyon has been moved, created, annihilated or 
fused when we apply the generalized Pauli operators (as we did in figure 10). This is done by analyzing the initial 
and the final state after the action of a bump operator and seeing which would be the energy after and before the 
process, as shown in figure 11. Therefore, we have: 



af' := (<T|)-1P|(,) + + a^P^^,^{a^)-' P^^^. (43) 



Here the upper-indices of operators aj are related to the energy cost of the process: 

• a^p"^ creates a pair of anyons of 2:- type and a^p annihilates it. The energy cost is AE" = 3. 

• and are related to the process of fusion or separation, respectively, of anyons as in figure 10.1 and 
also to the process of creation (and annihilation) of a pair of anyons tied to a previous pair. The energy cost is 
AE = l 

• moves anyons and also it can invert the orientation of a pair of anyons (as in figure 11.2). There is no energy 
cost in these processes. 

For the plaquette operators Bp we proceed in the same way obtaining a similar result. The corresponding local 
projectors that we denote as Rj are built analogously just by changing As for Bp^ where p and p' are the adjacent 
plaquettes to the qutrit j. Then the operators which describe the analogous process for x-type anyons are: 



(^i + ^lP+{2) + ^0(l)(^j )~^-^0(2)' 



(44) 



-(2)- 



Some of these operators are associated to more than one projector, unhke for qubits. That is because for 3-level 
systems, the possibihties for different excitations scenarios have grown significantly. 

As we have seen in the previous section, these operators arise naturally as the Fourier transform of the interaction 
Hamiltonian when a thermal bath is weakly coupled with our system, 
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5„e-«^°^^=^S„He--*. (45) 



In this case the interaction Hamiltonian will be of the form: 

= E 5„ ® /„ = ^ <t| ® /; + (<t|)-i ® [0 + a] ® /; + (<tJ)-i ® (/;)t, (46) 

and it is quite important to remark that there are only 3 Bohr frequencies this time, uj = 0, ±|, ±3. 

We can check that the dynamical operators obtained are indeed compatible with this interaction potential as 
'^oc = X^q: OUT casc, it is trivial to check: 



n 



with n = 0, 1, 2, using equations (44) and (45). 

Moreover, [iJ, a^] oc a^, based on the fact that H^^^ is made of stabilizers which at most introduces a phase when 
they are applied to states a* Thus, As{Bp)a'^ |^) oc a* |^) and a'^As{Bp) |(/)) oc a* |0) therefore [i^, a*] oc a*, Va* a 
dynamical operator of our system. With this proviso, the Davies generator turns out to be given by: 

^=Q{X)=iS{X)+C{X), (47) 

with 

5{X) = [H^y^,X] = l[H,^^ + Hl^,X], 
jC{X) = C'{X)+£^{X), 

C%X) = Y: \R{S){{-bf^'bf^X - Xjf'lf + 2hf^Xhf) + e-'^-hfhf^'x - Xbfhf^' + 

3 

+ 2bf^Xbf^'}) + \R{3/2){i-bfhfx - Xbfhf + 2bf'Xbf) + e-l^{-bf¥f^' X - 
- Xbfbf^ + 26f X6f t)} - ii?(0)[6°, [blX]], 



C^{X) = \R{m-4^'4^X - Xa^-'^a^ + 2a^'Xa^'>) + {-a^'> a^^ X - Xaf^ 



2--\-/L\ -J -J -- J -J ■ — J J / ■ ~ \ "3 "3 " 3 3 

3 

,(1) V.(l)tx^ , i P/Q /o^ „(2)t„(2) X. _ v.(2)t„(2) , o.(2)t J. .-|^r_.(2).(2)t ■ 



+ 2a'^'Xa'^>'}) + ^-Rm{{-a'f"aY>X - Xa>f"a'^' + 2a'^" Xa^') + e-^-P{-a'^' a^" X - 
- Xafaf^ + 2afXaf^)} - Ji?(0)[a°, [a°, X\]. (48) 



D. Topological Order 

Similarly to the case of qubits, we will study the evolution of the expectation value (GS| Xc |GS), where Xc is the 
tensor product of generalized Pauli operators {d = 3) along a non-contractible loop, and |GS) denotes a certain 
ground state in the stabilizer subspace; namely a superposition of the degenerate states in the ground state manifold 
of H^y\ 

In the weak coupling limit, the master equation that describes the dynamics of this quantity is: 

'^^=i[H^y^,xm+axm. m 

In order to simplify the calculation we remove the free evolution part of the equation 

XS) = e--^"''''X,{t)e^"'''' =^ ^ = C[XM, (50) 
being both the dissipator C and the mean value (GS| Xc |GS) invariant under this transformation. 
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E. Short time regime 

In the short time regime, we can approximate Xc{t) ^ {l-\-tC)Xc; here we denote Xc := Xc{0). Thus, the evolution 
of (GS|Xe(t)|GS) is 

{Xc{t)) ^ (GS| X, ICS) + t (GS| C{X,) ICS) . (51) 

We need to calculate (GS| C{Xc) |GS), with C{Xc) = C^{Xc) + C^{Xc). This calculation is made in Appendix A, 
obtaining 

(GS| |GS) = -|i?(A)e-^^|c| (GS| |GS) . (52) 

Hence, we can define F := yi?(A)e~^^|c| as the initial decay rate of the system. For qutrits A = 3 while for qubits 
(see Eq. (25)) we have obtained an analogous expression but with A = 4 instead. 

This result can be generalized for the case of qudits with arbitrary d. We have already seen, that at short times, 
only the creation of anyons contributes to the decay of topological order. The free diffusion of anyons and the fusion 
processes among them will not appear as they are second order processes in time. However, as we increase d there 
are more types of anyons with different energies. Moreover, a pair of anyons should always be compatible with the 
conditions fls = 1 Hp = 1- That means that the possible types of anyons with different energies are of 
the form uj'^ — uj^~^ with n = 1, [|J, and respective energies A^ = 2(1 — cos ^^). Note that n = 1 refers to the 
lowest energy pair of anyons, i.e. the energy gap of the Hamiltonian. Thus, the initial decay rate has to be the sum 
of all these contributions: 

L-J 

r. = E%l^l^(^n)e-^^^. (53) 

n=l ^ 

It is important to point out that in the case of qudits, an analogous expression for the interaction with the environment 
to (13) involves 5'^ = cr^, (cr^)^, (cr^)^~^, a^, (cr^)^, (<j^)^~^. Ah non trivial powers of and are included to 
allow for excitations of physical qudits from one level to another, at first order in time. 

Using Eq. (53) it will be possible to establish a crossover temperature Tc as the limit for which the initial decay 
rate F will be larger for qubits than for qudits. For the sake of comparison we take i^(A^) the same for qubits and 
qudits. This is reasonable since A^ are of the same order, and R{An) are the Fourier transforms of the bath coupling 
that induces the excitations on the physical qudits. Thus, we set up the condition Fc^(Tc) := F2 (Tc). Using Eq. (53) 
we arrive at the following expression: 

L-J L-J 
4=f]A„e-(^"-^)''>f]A„, (54) 

n=l n=l 

as A^ < 4 for d > 2, Vn. Therefore, this equation only has a solution for such values of d satisfying Xl^Li A^ < 4. 
But this is only true for d = 3. Thus, there exits only such a Tc for qutrits. For other values of d, the initial decay 
rate for qudits will always be larger than for qubits. This happens as J2n ^ri increases almost linearly with d, and 
= 3 is the only case when this quantity is smaller than 4, i.e., the gap in the case of qubits. Let us now compute Tc 
for qutrits: 



3Eoe-^^'^^ = 4Eoe-^^'^^ , (55) 
with Eq the natural energy unit of the system. This leads to the following crossover temperature, 

Te=^. (56) 
kB In I 

The meaning of this temperature is the following. Above this temperature Tc, the initial decay rate for qutrits is 
smaller than for qubits, something that makes qutrits better in this comparison. For Eq ^ lOOkHz used in the proposal 
of a Rydberg quantum simulator [64] for the operators of the 2-D Toric Code, we obtain an estimate of Tc ~ 20fiK. 
In addition, it could be computed a Tc comparing systems with d odd and {d— 1) even. There is always a temperature 
above which the system of qudits with d odd has a smaller initial decay rate than the previous {d — 1) even. 
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It is also important to point out that F is only the initial decay rate. It is possible that the dynamics of anyons, 
with free diffusion etc., play an important role in the loss of topological order. Beyond short times, our conjecture is 
that the new processes that appears in the case of qutrits, i.e. fusion of anyons which end tied up, will be an obstacle 
for the free diffusion of anyons. This would represent an improvement for the stability of the generalized toric code 
in some intermediate time regime for this is the cause of the loss of topological order in the system. 

F. Long-Time Regime 

Now we want to study the master equation (49) in the opposite time regime. We are interested in the fate of the 
non-local order parameter we are using to describe the topological order in a system of qudits in a generalized toric 
code. We conjecture that the final state will be given by a thermal Gibbs state. To show that our observable for the 
order parameter {Xc) approaches to the expectation value of Xc in the Gibbs state for times long enough, we resort 
again to the condition (7). In the generalized case, it reads as follows 

{a„ a^, . . . , (Jt-\(J,, (jI,..., G^-^y = CI, for any d. (57) 

This is due to the fact that if some generic operator, say A, commutes with every element of the set 
{cr^, cr^, . . . , a^~^ ^az^cTx^ • • • , crf~"^}, so does with every element of the d-Pauli group. This follows from the Jacobi 
identity and the fact that dzcrx = ^o-xCTz- Therefore given the irreducibility of the computational representation the 
(i— Pauli group (the technical details of this proof are given in Appendix B) the condition (57) holds. 
With this result, we may obtain the behaviour in the long time regime 

{X,{t ^ <^)) = MX,p{t ^ oo)) = I ^e-^^^ {iPi\X, IVi) = 0, (58) 

i 

which implies that the topological order is also destroyed for qudits in the generalized toric code when times of 
interaction with a thermal bath are long enough. 

Now, let us summarize and combine the results for both time regimes, i.e., short and long time behaviours. We 
have proved that at short times the global order parameter we are considering behaves as: 

{X,{t))p=e-^\X,{0)), (59) 

with r = ^R{A)e~^^\c\ and A = 3 for qutrits. We have also shown that there exits a crossover temperature Tc 
above which, the initial decay rate for qutrits is smaller than for qubits. Furthermore, we have shown this event only 
occurs in the case of qutrits, as for other values of the initial decay rate is always larger than for qubits. On the 
other hand, far from this initial short-time regime, the topological order of the system decays to zero for times long 
enough. 

IV. CONCLUSIONS 

We have introduced the basic concepts of 2-D Kitaev Model for qubits as well as a generalization of the code for 
qudits, i.e. d— level systems with the main purpose of studying its decoherence properties due to thermal effects. To 
this end, we have coupled these systems to thermal baths in order to study the thermal stability within a quantum 
open systems' formalism, namely Davies' theory. 

The generalization of the toric code leads to new physics. Indeed, we have particularized for the case of qutrits 
and obtained very interesting results. First of all, new abelian anyons have arisen with novel braiding properties, i.e. 
new statistics by exchange of particles. For instance, let us move a pair of anyons around another pair who stays 
still. We would pick up a different phase, letting the first pair still and moving the other one around. Furthermore, 
new energy processes appear which are forbidden for qubits, being d = 3 the first non-trivial system where these new 
processes can be observed. Moreover, we present a master equation that describes the dynamics of any observable of 
the system coupled to a thermal bath, giving a complete description of the problem. 

We have proposed a new way to study thermal stability regarding the loss of topological order in the system. At 
short times, the system starts loosing its order with a certain decay rate that we are able to compute explicitly. We 
have checked that the system relaxes to the thermal state for any value of d, as it was expected. However, we have 
proved that above a certain crossover temperature, the initial decay rate for qutrits is smaller than the one from the 
original case for qubits. Surprisingly, this behaviour only happens with qutrits and not with other qudits with d > 3. 
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It would be very interesting to be able to generalize further this study to other topological codes [65-71] coupled 
to thermal baths by deriving appropriate master equations for them. Other challenges in this direction are to study 
thermal effects with non-abelian topological codes [72-78], higher dimensional codes [12, 79-89] and systems with 
topological order based on two-body interactions [90-93] , instead of many-body interactions in the Hamiltonian. This 
would facilitate the physical simulation of these topological quantum models [64, 94-100]. 
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Appendix A: Evolution of the Order Parameter for Qutrits 

In order to compute {GS\ C{Xc) |GS) (with C{Xc) = C^{Xc) + C^{Xc)), we need the expression of the system 
operators that appear in Eq. (48) which were defined previously in Eq. (44) and (45) . These operators are expressed 
in terms of some orthogonal projectors whose definition is given in Eq. (42). However, there are only two projectors 
which are relevant here, namely 

|GS) = |GS) and |GS) = |GS) , (Al) 

as the rest of them vanish when acting on the ground state. Remember that are the projectors associated with 
the stabilizers As and with stabilizers Bp. Moreover we have 

bf^\GS)=0, bf\GS)=0, &f^^|GS) = 0, bf^\GS)=0. (A2) 
Thus, after doing some simplifications on Eq. (48): 

(GS| |GS) = ^e-^^ ^ (GS| {2bfxj;^^ - bfbf^^X, - xj^hf^) |GS) = 

3 

= 2\c\ (GS| (aj + {(7^)-')X, |GS) = 0, (A3) 

as Xc |GS) oc |GS) but crj |GS) is orthogonal to |GS), and we have used the fact that [P^ q,Xc] = for every j, as 
these projectors are only functions of vertex operators. This is not true for i?^ q J ^ i-^- J belongs to the path 
where Xc is acting on. In that case, since a^a^{a^)~^ = Co^crj , we obtain a^Xc{(J^)~^ = ooXc for the string operator. 
In addition, by making use of 

a5^^|GS) = 0, af |GS)=0, af ^ |GS) = 0, |GS) = 0, (A4) 
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the result for (GS| C^{Xc) |GS) turns out to be 

(GS| |GS) = = ^e-^^ (GS| i2af^xjp' - af\f^' - X^af\f^') |GS) = 

3 

= ^6-3/5 (GS| (^1 + (^|)-^)^cK + [a])-') |GS) - (GS| P^+X, |GS) - 

j 

- I (GS| Pi+(a| + (a|)-i)Pi+X, |GS) - (GS| X^P^^ |GS) - 

- (GS| XePi+(<7| + K)-i)Pi+ |GS) = 

= ^6-3/5 ^ 5,Vc((GS| (2 + a| + |GS) - (GS| X, |GS) - 

j 

- I (GS| (<7| + (<7|)-i)Xe |GS) - (GS| |GS) - I (GS| (<7| + (<7|)-i)Xe |GS)) + 
+ Sjec{{GS\ {a] + + |GS) - 2 (GS| X^ |GS) - 

- (GS| {a] + (<t;)-1)X, |GS)) = -^P(3)e-3^|c| (GS| X, |GS) = 

= -|i?(A)e-^''|c|(GS|Xe|GS), 
where |c| is the number of points in the path c. 



Appendix B: Irreducibility of the Computational Representation of the Pauli Group 

The (i-Pauh group is generated by products of and Oy such that cr^ = erf = 1 and Gz<^x = (^cTx^^z where uj is 
a primitive d-ioot of the unity. Its order is d^, which is a direct consequence that any element of the group can be 
written as cj^cr^cr^ for some n, m and k. 

We take the representation of the c?— Pauli group when acting on the computational basis: 

c^xItt) = 1^ + 1) mod. (i, (Bl) 
a,\n) = uj^\n), (B2) 

and we want to show this representation is irreducible. We proceed by computing the character x of every of its 
elements, which is given by the trace of the matrices. Using the computational basis when taking the trace, from 
the above relations, xi^T) = for m G {1, . . . , — 1}. Similarly xi^T) = for m G {1, . . . , — 1} as the sum of 
the roots of the unity vanishes. On the other hand, because azCJx = (^cTx^^z and the cyclic property of the trace, we 
conclude that the character of every element of the form cr^cr^ is zero for any representation. The rest of the terms 
are proportional to the identity cj^l, and so x(co'^l) = u^d. 

The irreducibility criterium asserts [9, 101] that a representation of a group G is irreducible if and only if the scalar 
product of characters is the identity, this is 

{X,X) = T^^YX*{9)X{9) = 1, (B3) 
1^1 sec 

where \G\ is the order of the group. For the computational representation of the c?— Pauli group we have 
thus, the representation is irreducible. 
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